The Surface Potential at the Air-water Interface Computed Using Density Functional 

Theory 
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An accurate prediction of the surface potential (cf>) at the air-water interface is critical to calcu- 
lating ion hydration free energies and electrochemical half-cell potentials. Using Density Functional 
Theory (DFT), model interfacial configurations, and a theoretical definition of cj>, we report a value 
of +3.63 volt at 0.92 g/cc water density. A maximally localized Wannier function analysis confirms 
that <f> is dominated by molecular quadrupole (or "spherical second moment") contributions. We find 
that the predicted surface potential depends on computational details, and conclude that standard 
DFT codes and the existing definition of <f> may not yield surface potentials directly comparable to 
existing experiments. 



"Surface potentials" (or "potentials of the phase" ) are 
invoked when a charged particle moves from one phase 
to another through their mutual interface. One of the 
simplest and most important examples pertains to the 
surface between air and water ("Air," "vacuum," 
and "vapor" will be used interchangeably in this work.) 
The electrostatic potential difference (4>) between an ion 
in vacuum at infinity and that ion in pure liquid wa- 
ter contributes a term q<f> to the absolute hydration free 
energy (AGhyd), where q is the ionic charge. <j> thus 
contributes to all aqueous media electrochemical half-cell 
potentials ; 17 ^ 8 which consist of AGhyd plus the pertinent 
ionization potentials. Recent Density Functional The- 
ory (DFT) calculations have yielded <j> estimates for the 
air- water interface 1 ^— substantially larger in magnitude 
than classical force field predictions^— and experimental 
values^ 1 - 2 -"— suggesting that the methodology and physi- 
cal basis of such calculations should be re-examined. 

We first focus on molecular simulation perspectives^— 
deferring the question of what is measured in experi- 
ments to the concluding discussions. For this purpose, 
we consider a charge-neutral slab of salt-free liquid wa- 
ter. Using this kind of simulation cell, force field-based 
simulations (i.e., those not based on electronic structure 
calculations) have reported that cj> depends on the wa- 
ter model usedfSr— but that the widely applied SPC/E 
model 24 yields ^=-0.55 volt, 8 similar to related three- 
point point-charge water models, including polarizable 
ones 3 

DFT takes into account electronic structure and molec- 
ular polarizabilites, and might be considered an improve- 
ment over non-polarizable water models. As will be 
shown, however, care must be exercised when interpret- 
ing the results based on the existing theoretical defi- 
nition of the surface potential (see below) £r— A DFT 
work that extrapolates the highest occupied molecular 
orbital in gas phase ion/water-clusters to infinite cluster 
size has estimated that <j> ~ 4 voh— (in our notation). 
Using the VASP codej 2 ^ the Perdew-Burke-Ernzcrhof 
(PBE) functional,— a bulk liquid simulation cell, and ap- 
proximations for surface dipoles, we have deduced that 
4> = +4.05 volt<21 By considering pK a and proton hy- 



dration free energies in DFT simulations, Cheng et al£l 
have postulated that (j> — +3.5 volt. With explicit air- 
water interfaces, another DFT calculation has predicted 
</> = +3.1 volt M In this work, we perform a DFT/PBE 
<f> calculation with air-water interfaces to compare with 
previous predictions. We also use maximally localized 
Wannier functions 2 ^ to analyze the results in terms of 
molecular contributions, which allows a detailed compar- 
ison of DFT and force-field based work.i^r— Our work 
highlights the dependence of (f> on computational details. 

The theoretical <fi is given by the difference in the aver- 
age plateau values between the liquid and vacuum regions 
in <j)(z), where 



dxdy V(x, y, z)j A 



(1) 



V(r) is the calculated electrostatic potential at point r, 
A is the lateral area of the simulation cell, z is perpen- 
dicular to the interfaces, and the angular brackets de- 
note statistical averaging. The analytic expressions de- 
rived for surface potentials, and formulas for dealing with 
long-range electrostatics within periodic boundary con- 
ditions in general, are based on purely coulombic, 1/r 
potentials^i 2 ^ Thus 



V(v) = J drV e (r')/|r-r'| + EVI 



R, - r| 



(2) 



Here p e (r) is the valence electron density, is the po- 
sition of pseudo- nuclei i, and Zi is the pseudo- nuclear 
charge in the pseudopotentials (PP), with Zo=+6|e| and 
Zh= +|e|. By default, VASP instead outputs the nega- 
tive of Eq. [2] after replacing the last term with the entire 
local PP J7(Rj,r). The short range, non-coulombic con- 
tribution to U(Ri, r) should be removed when computing 
cj). This ambiguity arises from the use of PP's, and is ab- 
sent in all-electron calculations. 

100 water configurations are taken from a one- 
nanosecond water-vapor interface molecular dynamics 
trajectory generated with 128 SPC/E water molecules 
in a 50.0x12.5x12.5 A 3 simulation cell (see the Sup- 
porting Information, SI, for more details). 4>{z) is com- 
puted from these configurations using DFT /PBE and the 
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FIG. 1: (a) Electrostatic potential (4>(z)) computed perpen- 
dicular to water-vacuum interface direction z, averaged over 
the lateral (a;- and y-) directions. The red and green lines rep- 
resent predictions using purely coulomb interactions (Eq. [2]) 
and U(Ri, r) (see text). J fi{z) = over the simulation cell, 
(b) Quadrupole (red) and dipole (green) contributions to <fi{z). 
The latter has its finite electric field removed and is then ex- 
panded ten-fold; the quadrupole component contains more 
spatial noise than <j>(z). 



Eq. [2] definition. Figure [Th shows that <j> averages to 
+3.63±0.04 volt when referenced to vacuum. If we had 
used the VASP default and included the short range part 
of the PP, (j> would become +2.8 volt instead. Clearly, 
the details of the algorithm strongly affect the results. 
The vacuum region exhibits a small, 0.005 volt/A elec- 
tric field, indicating that a small but finite net average 
dipole moment persists in the water slab due to insuffi- 
cient statistical sampling. 

<fi can be rigorously decomposed into quadrupole and 
dipole components ^UL2& In water, they can be written 

(fi = (fi q + (fid 

= ^g(^watcr) ~ 4>q{z a ir) + </>d (Zwatcr ) ~ ^d^eix)] 



-(2tt/3A) f dr<5(z-i? m ,,)$> m (r)(r-R, 
Uk/A /dr^p m (r)(r z - i? m , z )9(z - i? m , 



Here p m (v) = pf{v) + Y. l ^r<K r - R i) is the charge den- 
sity of molecule to with all electrons and nuclei i residing 
on rn, R m is the molecular center (oxygen atom in the 
case of water), and 9 is the Heaviside function. Calculat- 
ing (fid requires an interfacial geometry. As an interface 
was not present in Ref. 21, (fid was approximated using 
the SPC/E water value (+0.21 volt) 4 

In contrast, <fi q can be computed in a bulk liquid simu- 



lation cell, with electron density demarcated into molec- 
ular contributions using maximally localized Wannier 
functions. 2 ^ Predicted to be +3.85 volt at 1.0 g/cc 
water density,— this VASP/PBE cfi q is lar ger m magni- 
tude than and opposite in sign to that for SPC/E water 
{4> q = —0.76 volt) because of differences in the charge 
distributions (Fig. [2]) . Using the O atom as the molecu- 
lar center, only the partial positive charges on the H-sites 
of SPC/E water contribute to Eq. yielding a negative- 
definite (fiq. For VASP/PBE, Eq.U]is instead dominated 
by the valence electron cloud surrounding the O-nuclei 
and <p q changes sign. 

To further analyze the results, we also use maximally 
localized Wannier functions to decompose the interfa- 
cial 4>{z) (Fig. la) into quadrupole and dipole contri- 
butions (Fig. lb). In the bulk liquid region, defined as 
5 A < z < 20 A, (fi q averages to +3.50+0.01 volt. This 
differs from the +3.85 volt derived indirectly (see the 
SI) because unlike bulk SPC/E water calculations j2i our 
small lateral simulation cell dimensions dictate a small 
Lennard- Jones cut-off distance for SPC/E water, which 
reduces the density of the bulk liquid region to 0.92 g/cc. 
Equation 3] implies that (fi q is proportional to the liquid 
water densityjii Consistent with this formula, the 8% 
reduction in p wa tcr and the 9.6 % decrease in <fi q com- 
pared to the indirect calculation conducted at 1.0 g/cc 
water density indeed track each other. We have not used 
PBE-based molecular dynamics to generate water slabs 
partly because PBE exhibits water densities that deviate 
even more strongly from experiments^ The +3.1 volt 
DFT value reported previously 23 likely also reflects the 
low DFT water density present in that work. 

The small cell also affects (fid- After removing 
the finite average electric field in the 100 configura- 
tions selected, the SPC/E (fid amounts to +0.009 volt, 
strongly reduced from the ^=+0.21 volt computed 
in a larger box£ DFT/PBE applied to these SPC/E 
interfacial configurations tracks SPC/E results, giving 
^=+0.012+0.008 volt (Fig. lb). (fi d and (fi q thus add to 
+3.52 volt, which is consistent with the +3.63 volt ob- 
tained directly using Eq. 1. The small discrepancy may 
arise from the fact that the interfacial simulation cell has 
prevented the extrapolation of Wannier estimates of <fi q 
(8) infinite box size as was done in Ref. [2(| 
A The term "quadrupole moment" used to describe <fi q 
is taken from the liquid state literature jldiLs 
"radial" and "spherical second moment" 20 ' 29 have been 
used elsewhere. Equation [4] contains the trace of 
quadrupole tensor, typically set to zero in multiple ex- 
pansions and does not contribute to real-space electro- 
static interactions^ As an example, the atoms in a neon 
solid are largely spherically symmetric and exhibit only 
small multipole (including off-diagonal quadrupole) mo- 
ments, but at the equilibrium density of 1.444 g/cc at its 
melting point, our DFT/PBE/PP (fi q estimate for solid 
neon still amounts to 3.6 volt. 

This neon example emphasizes that the theoretical <fi 
in Eq. [TJ computed using the PBE functional, is not a 
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FIG. 2: Charge distributions, (a) SPC/E water model; (b) 
projection of DFT/PBE water, using pseudopotentials. O, 
H, and electron density exceeding 0.05 |e|/A 3 are depicted in 
red, black, and orange respectively. 



physical quantity. Indeed, the PP used omits oxygen 
core Is electron contributions (see the SI); otherwise the 
magnitude of <j) would be even larger (Eq. This 
theoretical <f> is however critical for DFT-based molec- 
ular dynamics calculations of ion hydration free energies 
(AGhyd)^ 1 ^ The reason, familiar in the classical force 
field literature jiJtiSiil is that AGhyd at infinite dilution 
are most conveniently computed as intrinsic ion hydra- 
tion free energies calculated in bulk- water boundary con- 
dition simulation cells using Ewald summation, with the 
non-ion-specific pure water "surface potential" contribu- 
tion q<f> added during post-processing. Ewald tcchnques 
arbitarily set the average electrostatic potential to zero 
over the bulk- water simulation cell. The "true" (but com- 
putation protocol-specific) average electrostatic potential 
is restored by referencing the liquid region plateau <f)(z) 
value to vacuum. This of course involves a rigid shift 
equal to 4>, which must be obtained using the same defini- 
tion and pseudopotentials applied in intrinsic ion hydra- 
tion calculations. In DFT calculations, <fi and the intrin- 
sic ion hydration free energy (the latter through Ewald 
summation conventions) both contain large, equal but 
opposite contributions from the water atomic core regions 
(Fig. [5J))2£ even though ions do not penetrate into water 
nuclei. cf> must be on the order of +4 volt to yield AGhyd 
comparable to experimental data.= If the VASP/PBE <fi 
were -0.55 volt like in SPC/E water, hydration of Cl~ 
would have been endothermic and unphysicalf 2 - Work 
function calculations in metals take advantage of a simi- 
lar cancellation of ambiguities^ 

We have so far side-stepped the issue of comparison 
with experiments. To the extent that the theoretical <p is 
mainly of interest for computing the absolute hydration 
free energies of ions at infinite dilution, it can be treated 
as a method- or force field-dependent entity. Indeed, it 
has been argued that the Galvani potential difference be- 
tween two phases is extremely difficult to measure, al- 
though possible in principled <j) also depends non-trivially 
on the salt present at the surfaced To our knowledge, in 
the experimental literatures^— there has not been a pre- 
cise definition of <j) in terms of microscopic (i.e., electronic 



and ionic) properties^ To the extent that this quantity 
has been indirectly measured at the air-water interface, 
early reported values strongly varied in magnitude and 
signp^ Several post-1970 experimental values are in bet- 
ter agreement with each other ; 12 ' 14 yielding +0.025 to 
0.16 volt values (see Ref. [TH, Table 1). They are also 
in reasonable agreement with the dipole contribution 4>d 
of the SPC/E water models which raises the intriguing 
point that these measurements may predominantly re- 
flect 4>d- However, these experimental values cannot be 
directly used to help calculate DFT-based ion hydration 
free energies for reasons discussed above. 

It may become possible to establish unambiguous ex- 
perimental surface potentials in the future. Here we con- 
fine ourselves to the observation that, just as the the- 
oretical <p depends on the method used, the measured 
(j) may be sensitive to experimental details. (1) Surface 
sensitive spectroscopic measurements may yield values 
for <j) that most closely match all-electron (frozen oxygen 
Is electrons or otherwise) DFT results. They amount 
to using as probes test particles that are point charges 
and do not exhibit many-Fermion effects. Proposed mea- 
surement of <j) using electron reflectivity^ with 1-10 keV 
beams should also probe the nuclear region, although 
many-Fermion effects may arise. (2) If ions are used 
as experimental probes, e.g., by considering ion hydra- 
tion free energies, the solvent atomic nuclear regions are 
not sampled. These regions are mainly responsible for 
the large and positive cj) in DFT calculations. Therefore 
a cf> value different from case (1) above should emerge. 
It is conceivable that cj) predicted using classical force 
fields, which do not contain distributed charges in atomic 
nuclear regions, may be more appropriate here. (3) In 
electrochemical measurements electrons are added 
to/removed from electrodes while ions enter/depart elec- 
tric double layers. Hence a mixture of electrons and ions 
are implicitly used as probes. In light of the present work, 
the results should contain <f> contributions from atomic 
nuclear regions of electrodes but not from water. 

In conclusion, using SPC/E H2O model-derived 
interfacial configurations, we have shown that our 
DFT/PBE pseudopotential calculation directly yields a 
+3.63±0.04 volt surface potential when we use the the- 
oretical definition of (f> and a water slab with interior 
density of 0.92 g/cc. Maximally localized Wannier func- 
tion analysis confirms that this value is dominated by 
the density-dependent molecular quadrupole (or "spher- 
ical second moment") contribution <f> q , consistent with 
previous calculations ! 20 ' 21 This suggests that the most 
robust way to estimate the dominant <p q term in DFT 
calculations is to use a bulk liquid water simulation cell 
at 1.0 g/cc density. The DFT surface dipole component 
of <j) is not accurately determined due to finite size ef- 
fects, but it tracks the SPC/E value for this simulation 
cell. This theoretical <f> is not a physical quantity, in the 
sense that it is not what is measured in existing electro- 
chemical experiments*^ - — But it is critical for calculatin 
ion hydration free energies and modeling electrochemical 



half cell reactions. In general, <j) computed in DFT cal- 
culations are sensitive to simulation details and may not 
be directly compared to experiments. 
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